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ABSTRACT 

Numerical  simulations  are  performed  using  the  Fluent  Computational  Fluid  Dynamics  code 
and  the  various  turbulence  models  to  simulate  the  horseshoe  vortex  formed  in  a  typical  wing- 
junction  turbulent  flow  experiment.  Simulations  were  conducted  using  the  renormalizable  k-s 
model,  the  Reynolds  Stress  Model,  the  V2F  model,  the  Spalart-Allmaras  model  and  the  k-co 
model.  The  calculated  results  were  compared  with  experimental  results  obtained  from  an 
extensive  database  available  on  the  internet  [3].  The  realizable  k-s  model  was  noticeably  less 
accurate  than  all  other  models  in  simulating  the  mean  velocity  components,  while  the 
remaining  models  all  displayed  similar  levels  of  accuracy.  None  of  the  models  were  able  to 
accurately  simulate  the  correct  behaviour  of  the  mean  kinetic  energy  as  a  function  of  position. 
The  V2F  model  however  came  closest  to  predicting  the  correct  behaviour,  and  offers  the  best 
combination  of  computational  accuracy,  computational  efficiency,  and  ease  of  use. 
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Simulation  of  a  Wing-Body  Junction  Experiment 
using  the  Fluent  Code 


Executive  Summary 

The  horseshoe  vortex  formed  at  the  sail/ hull  junction  on  a  submarine  has  a  number  of 
characteristics  which  are  undesirable  for  efficient  submarine  performance.  These 
include  the  generation  of  unwanted  noise,  increased  local  surface  shear  stress  and 
drag,  as  well  as  incoherent  fluid  loading  on  the  propellers  which  can  lead  to 
narrowband  sound  generation.  An  increased  understanding  of  the  fluid  behaviour  in 
these  types  of  flows  would  lead  to  improvements  in  submarine  performance  through 
increased  stealth,  reduced  vulnerability  and  improved  hydrodynamic  performance. 

A  detailed  experimental  study  of  the  horseshoe  vortex  formed  in  the  vicinity  of  the 
nose  of  a  modified  NACA  0020  airfoil  in  a  fully  developed  turbulent  boundary  layer 
has  recently  been  reported  in  the  literature  by  Devenport  and  Simpson.  The  data 
obtained  in  this  experiment  are  of  a  particularly  high  quality  and  are  readily  available 
in  electronic  format  on  the  internet.  This  report  documents  results  obtained  using  the 
Fluent  Computational  Fluid  Dynamics  code  to  simulate  the  horseshoe  vortex  formed  in 
this  experiment.  It  forms  part  of  an  overall  study  aimed  at  understanding  and 
characterizing  the  vortex  near  the  sail  and  at  downstream  locations  so  that  designs  may 
be  implemented  to  minimize  the  adverse  effects  of  the  vortex  structure  on  submarine 
performance. 

Simulations  were  conducted  using  five  different  turbulence  models  contained  within 
the  code  to  ascertain  the  effectiveness  of  each  model.  The  calculated  results  were 
compared  with  experimental  results  obtained  from  the  database  available  on  the 
internet.  The  simulations  showed  that  the  realizable  k-s  model  was  noticeably  less 
accurate  than  all  other  models  in  simulating  the  mean  velocity  components.  The 
remaining  models,  which  consisted  of  the  Reynolds  Stress  Model,  V2F  model,  Spalart- 
Allmaras  model  and  the  k-co  model,  all  displayed  similar  levels  of  accuracy.  None  were 
able  to  accurately  simulate  the  correct  behaviour  of  the  mean  kinetic  energy  as  a 
function  of  position  close  to  boundaries,  although  all  were  able  to  simulate  the  time- 
averaged  velocity  components  at  most  locations  within  the  flow  to  within  an 
acceptable  accuracy.  The  V2F  model  came  closest  to  predicting  the  correct  behaviour 
overall,  and  offers  the  best  combination  of  computational  accuracy,  computational 
efficiency,  and  ease  of  use  for  further  studies  of  this  type  of  turbulent  flow. 
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1.  Introduction 


The  Hydrodynamics  Group  within  Maritime  Platforms  Division  (MPD)  performs  computational 
fluid  dynamics  (CFD)  studies  on  underwater  flow  problems  relevant  to  submarine  design.  These 
studies  examine  and  aim  to  improve  the  performance  of  these  platforms  through  increased  stealth, 
reduced  vulnerability  and  improved  hydrodynamic  performance.  Problems  of  interest  involve  the 
simulation  of  flow  around  a  manoeuvring  submarine,  hydrofoil  studies,  hull/ propeller 
interactions,  viscous  bluff  body  modelling,  stores  release  modelling,  and  the  study  of  the 
horseshoe  vortex  formed  in  the  boundary  layer  at  the  hull/ sail  junction.  These  studies  are  carried 
out  using  a  variety  of  CFD  codes  including  the  finite  volume  codes  Fluent  [1]  and  CFX  [2]. 

The  horseshoe  vortex  formed  at  the  sail/ hull  junction  on  a  submarine  has  a  number  of 
characteristics  which  are  undesirable  for  efficient  submarine  performance.  The  vortex  is  subject  to 
large-scale  low-frequency  bistable  unsteadiness  which  appears  to  be  generated  in  the  nose  region 
and  which  leads  to  the  generation  of  unwanted  noise.  The  vortex  entrains  high-momentum  free- 
stream  fluid  into  the  corner  between  the  hull  and  sail  and  greatly  increases  the  local  surface  shear 
stresses  and  drag.  The  trailing  legs  of  the  vortex  persist  very  far  downstream  of  the  sail  and  can 
maintain  coherence  as  they  are  ingested  into  the  propellers.  This  provides  a  source  of  unsteady 
loading,  whose  adverse  side  effect  may  be  narrowband  sound  generation. 

Devenport  and  Simpson  [4]  carried  out  a  detailed  study  of  a  horseshoe  vortex  in  the  vicinity  of  the 
nose  of  a  modified  NACA  0020  airfoil  for  a  fully  developed  turbulent  boundary  layer  and  were  the 
first  to  report  the  bistable  nature  of  the  flow.  Fleming  et  al.  [5],  Olcman  and  Simpson  [6]  and 
Devenport  et  al.  [7]  conducted  further  experiments  on  the  same  wing  geometry  to  further  elucidate 
details  on  the  formation  and  structure  of  the  vortex  near  the  appendage.  The  velocity  field 
downstream  of  the  appendage  has  not  been  studied  to  the  same  extent  however,  and  to  address 
this  shortfall  Anderson  and  Huyer  [8]  recently  used  particle  image  velocimetry  (PIV)  to  measure 
the  full  three-dimensional  mean  and  turbulent  velocity  field  at  various  stations  downstream  of  a 
standard  NACA  0020  airfoil. 

This  report  documents  results  obtained  using  the  Fluent  code  to  simulate  the  horseshoe  vortex 
formed  in  the  Devenport-Simpson  experiment  [4].  It  forms  part  of  an  overall  study  aimed  at 
understanding  and  characterizing  the  vortex  near  the  sail  and  at  downstream  locations  so  that 
designs  may  be  implemented  to  minimize  the  adverse  effects  of  the  vortex  structure  on  submarine 
performance.  The  advantages  of  studying  this  particular  experiment  are  the  extensive  set  of  high- 
quality  measurements  obtained,  the  ready  availability  of  this  data  in  electronic  format  [3]  and  the 
detail  with  which  the  turbulent  boundary  layer  in  the  approach  flow  has  been  characterised. 
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2.  The  Devenport-Simpson  Experiment 


Devenport  and  Simpson  [4]  studied  the  behaviour  of  a  turbulent  boundary  layer  on  a  flat  plate  as 
it  encountered  the  nose  of  a  cylindrical  wing  mounted  normal  to  the  plate.  The  wing  cross-section 
consisted  of  a  3:2  elliptical  nose  (with  its  major  axis  aligned  with  the  chord)  and  a  NACA  0020  tail 
joined  at  the  maximum  thickness.  The  chord  length  was  30.5  cm,  maximum  thickness  7.17  cm  and 
the  height  was  22.9  cm.  The  wing  was  mounted  on  the  test  wall  of  an  open-circuit  type  wind 
tunnel  powered  by  a  centrifugal  blower.  The  test  section  was  91  cm  wide  and  25.5  cm  high.  The 
boundary  layer  on  the  test  wall  was  tripped  by  a  0.63  cm  high  step  prior  to  entering  the  test  section 
so  that  the  inlet  conditions  were  clean  and  well  defined.  The  shape  of  the  wind  tunnel  leading  into 
the  test  section  was  such  that  the  streamwise  pressure  gradient  was  effectively  zero,  thus 
producing  an  equilibrium  turbulent  boundary  layer  on  the  test  wall.  The  turbulence  intensity  was 

0.2%  at  the  reference  speed  (Uref)  of  27  m/s  and  the  value  of  the  Reynolds  number,  based  on 
momentum  thickness,  was  6700.  Figure  1.  shows  a  schematic  of  the  test  geometry. 

A  right-handed  coordinate  system  was  used  where  X  is  in  the  direction  of  the  free-stream  flow,  Y 
is  normal  to  the  flat  plate  and  Z  is  normal  to  the  wing  symmetry  plane.  The  origin  of  the 
coordinate  system  is  at  the  upstream  foot  of  the  wing  and  T  is  the  maximum  thickness  of  the  wing. 
The  corresponding  time-averaged  velocity  vectors  are  denoted  by  U,V,W  and  the  instantaneous 
velocities  by  u,v,w.  A  three-component  laser  Doppler  velocimeter  was  used  to  measure  all  non¬ 
zero  mean-velocity  and  Reynolds-stress  components  at  eleven  stations  in  the  plane  of  symmetry 
upstream  of  the  wing.  The  most  upstream  station  was  located  at  X/T  =  -0.86  and  the  closest  station 
to  the  wing  was  located  at  X/T  =  -0.05.  At  each  station  measurements  were  made  along  a  line  in 
the  Y  direction  covering  the  range  0.001  <  Y/T  <  1.0. 


X=1.3m 


Z=0.5m 


Figure  1.  Perspective  view  of  the  wing-body  junction. 
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The  adverse  pressure  gradient  in  the  streamwise  direction  imposed  by  the  presence  of  the  wing 
causes  the  approaching  boundary  layer  to  separate  and  form  a  horseshoe  vortex  which  wraps 
around  the  wing.  The  experiment  by  Devenport  and  Simpson  studied  both  the  time-dependent 
and  time-averaged  turbulence  structure  near  the  nose  of  the  wing-body  junction.  Their  time- 
averaged  velocity  measurements  showed  a  single  vortex,  roughly  elliptical  in  shape,  centred  near 
X/T  =  -0.2,  Y/T  =  0.05.  The  vortex  generated  an  intense  backflow  with  a  maximum  mean  velocity 

of  0.48Uref  and  then  decelerated  to  the  separation  point  at  X/T  =  -0.47.  From  the  most  upstream 
station  to  the  vicinity  of  the  separation  point  the  measurements  showed  a  flow  developing  in  a 
manner  similar  to  that  of  a  two-dimensional  turbulent  boundary  layer  in  an  adverse  pressure 
gradient.  In  the  vicinity  of  the  junction  vortex  however,  between  X/T  =  -0.3  and  -0.15,  the  profiles 
of  the  mean  and  turbulence  quantities  ceased  to  bear  much  resemblance  to  those  in  an  adverse 
pressure  gradient.  In  this  region  the  turbulence  stresses  were  much  greater  and  reached  values 
many  times  larger  than  those  normally  observed  in  turbulent  flows.  Examination  of  the  time- 
dependent  data  showed  that  these  large  stresses  were  associated  with  bimodal  (double-peaked) 
histograms  of  velocity  fluctuations  produced  by  a  bistable  velocity  variation.  Figure  2  shows  a 
schematic  of  the  mean-time  velocity  vectors  in  the  plane  of  symmetry  which  illustrate  the  overall 
nature  of  the  flow. 


Figure  2.  Mean-time  velocity  vectors  in  the  plane  of  symmetry 

Simpson  [9]  has  discussed  the  time-dependent  nature  of  the  flow  in  the  junction  region  in  detail 
and  summarises  much  of  the  experimental  work  on  such  systems  prior  to  the  detailed  experiments 
of  Devenport  and  Simpson.  Initially  one  large  horseshoe  vortex  exists  in  the  junction  region.  As 
this  vortex  is  stretched  and  weakened  a  secondary  separation  vortex  is  formed  downstream  of  the 
separation  and  sometimes  even  tertiary  vortices  can  be  formed.  These  vortices  appear  to  merge 
with  the  original  weakened  vortex  at  some  stage  and  temporarily  strengthen  it,  then  the  whole 
structure  becomes  unstable  and  is  replaced  by  a  new  large-scale  horseshoe  vortex  and  the  cycle 
repeats,  although  the  process  never  becomes  periodic. 

3.  Previous  Computational  Studies 

Apsley  and  Leschziner  [10]  reported  the  results  of  a  computational  study  of  the  predictive 
performance  of  twelve  turbulence  models,  within  three  distinctly  different  classes,  for  the  wing- 
body  junction  flow  experimentally  studied  by  Devenport  and  Simpson  [4]. The  turbulence  models 
consisted  of  five  linear  eddy  viscosity  models  (including  variants  of  both  the  k-z  and  k- oo  models). 
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three  non-linear  eddy  viscosity  models  (all  of  the  k-z  type),  and  four  Reynolds  Stress  Models.  The 
computations  were  performed  as  part  of  a  collaborative  turbulence-model-validation  study 
involving  six  partners:  UMIST,  Loughborough  University,  BAE  Systems,  DERA,  ARA  and  Rolls- 
Royce  Aeroengines.  Model  performance  was  judged  by  comparing  solutions  with  experimental 
data  for  pressure  fields  on  the  plane  wall  and  around  the  aerofoil;  for  velocity,  turbulence  energy, 
shear  stress  and  streamwise  normal  stress  in  the  upstream  symmetry  plane;  and  for  velocity, 
turbulence  energy  and  shear  stress  in  cross-flow  planes  downstream  of  the  aerofoil  leading  edge. 
The  main  finding  of  the  study  was  that  Reynolds  Stress  Models  offered  predictive  advantages  over 
the  other  models  examined,  especially  in  terms  of  the  far-field  structure  of  the  horseshoe  vortex, 
although  no  model  achieved  close  agreement  with  the  experimental  data  in  respect  of  both  mean 
flow  and  turbulence  quantities. 

Parneix  et  al.  [11]  used  the  V2F  turbulence  model  to  simulate  three  complex  three-dimensional 
turbulent  flows:  the  infinitely  swept  bump,  flow  in  a  U-bend,  and  the  appendage-body  junction 
flow  of  Devenport  and  Simpson  [4].  The  V2F  model  makes  use  of  the  standard  k-z  model  but 
extends  it  by  incorporating  near-wall  turbulence  anisotropy  and  non-local  pressure-strain  effects, 
while  retaining  a  linear  eddy  viscosity  assumption.  It  uses  fewer  equations  than  a  full  Reynolds 
Stress  Model  and  is  more  numerically  robust  than  these  models.  The  V2F  simulation  was  able  to 
accurately  calculate  the  separation  location  as  well  as  the  intensity  of  the  back-flow.  The  computed 
U  and  V  profiles  in  the  symmetry  plane  showed  excellent  agreement  with  experiment.  In  contrast 
with  all  the  k,z  and  RSM  models  surveyed  by  Apsley  and  Leschziner  [10],  the  V2F  model  was  also 
able  to  accurately  simulate  the  location,  shape  and  size  of  the  region  where  the  high  levels  of 
turbulence  existed,  and  the  calculated  values  for  k  were  only  about  30%  lower  than  the 
experimental  values. 

Rodi  et  al.  [12]  surveyed  the  results  of  computations  submitted  to  two  ERCOFTAC  Workshops  on 
Turbulent  Flow  Calculations  for  five  test  cases,  including  the  Devenport  -Simpson  experiment.  The 
calculations  were  made  using  a  wide  variety  of  turbulence  models  including  the  standard  k-z 
model  with  wall  functions,  various  low-Reynolds-number  versions  of  the  k-z  model,  the  RNG 
version  of  the  k-z  model,  various  non-linear  versions,  the  k-o o  model,  various  versions  of  Reynolds- 
stress  models  and  an  algebraic  stress  model.  The  general  conclusion  for  all  the  test  cases 
considered  was  that  the  more  complex  Reynolds-stress  models  and  intermediate  algebraic-stress 
and  non-linear  k-z  models  were  not  consistently  better  than  the  simpler  linear  eddy  viscosity 
models  for  the  practically  relevant  mean  quantities,  but  in  some  cases  they  were  found  to  give 
superior  predictions  of  the  details  of  the  turbulence  quantities. 

For  the  specific  case  of  the  Devenport-Simpson  experiment  five  calculations  used  versions  of  the  k- 
z  model  and  one  used  a  Reynolds-stress  model.  All  of  the  calculations  reproduced  reasonably  well 
the  pressure  distributions  on  the  ground  plate  and  the  wing  as  well  as  the  general  flow  pattern, 
although  for  most  models  the  size  of  the  vortex  was  predicted  to  be  too  small,  particularly  for  the 
RSM  model.  The  strength  of  the  horseshoe  vortex  as  it  swept  around  the  wing  was  also  under 
predicted  by  all  the  models.  Experimentally  the  highest  turbulent  kinetic  energy  occurs  in  the 
centre  of  the  vortex.  However  all  the  calculations,  except  those  using  the  RNG  version  of  the  k-z 
model,  predicted  the  maximum  value  of  k  to  be  very  close  to  the  corner,  with  excessive  production 
of  k  in  the  stagnation  region  in  front  of  the  leading  edge.  This  is  a  well  known  problem  with  the 
standard  k-z  model,  and  also  occurred  with  the  RSM  model.  Only  the  RNG  version  of  the  k-z 
model  did  not  have  this  problem  and  predicted  the  maximum  value  of  k  in  the  vortex  core, 
although  with  a  much  smaller  value  than  in  the  experiment. 

Chen  [13]  analysed  the  performance  of  two  turbulence  models:  a  two-layer  isotropic  eddy  viscosity 
model  and  a  near-wall  Reynolds-stress  model.  The  RSM  model  is  a  low  Reynolds  number 
extension  of  the  high  Reynolds  number  model  of  Speziale  et  al.  [14]  which  is  valid  all  the  way  to 
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solid  walls.  The  calculations  clearly  demonstrated  the  superior  performance  of  the  RSM  model 
over  the  simpler  isotropic  eddy  viscosity  model.  In  particular,  the  RSM  model  successfully 
resolved  the  detailed  turbulence  structures  associated  with  the  development  of  the  appendage 
boundary  layer,  the  formation  of  the  spiral  primary  horseshoe  vortex  and  the  corner  secondary 
horseshoe  vortex,  and  the  evolution  of  both  vortices  around  the  appendage/ flat-plate  junction. 

Fu  et  al.  [15]  studied  the  behaviour  of  a  standard  linear  eddy  viscosity  with  two  quadratic  non¬ 
linear  eddy  viscosity  models:  the  Gatski-Speziale  model  [16],  and  a  realizable  version  of  the 
Gatski-Speziale  model  due  to  Fu  et  al.  [17].  All  three  models  gave  close  results  for  the  mean  flow 
velocities,  but  significant  differences  appeared  in  the  Reynolds  stresses.  The  two  non-linear  eddy 
viscosity  models  provided  overall  better  performance  than  the  linear  version,  with  the 
renormalized  version  providing  marginally  better  performance.  The  three-dimensional  boundary 
layer  skewness  was  not  calculated  correctly  however  and  Fu  et  al.  concluded  that  three- 
dimensional  anisotropic  effects  could  not  be  modelled  using  a  quadratic  non-linear  eddy  viscosity 
model. 


4.  The  Fluent  Code 


Fluent  Version  6.1  is  a  CFD  computer  code  developed  and  marketed  by  Fluent  Inc.  [1].  The  code 
solves  the  equations  for  conservation  of  mass,  momentum,  energy  and  other  relevant  fluid 
variables  using  a  Finite  Volume  technique.  First  the  fluid  domain  is  divided  into  a  large  number  of 
discrete  control  volumes  (also  known  as  cells)  using  a  pre-processor  code  which  creates  a 
computational  mesh  on  which  the  equations  can  be  solved.  The  meshing  software  available  with 
Fluent  is  called  Gambit.  This  software  can  create  a  two-dimensional  (2D)  surface  mesh  using 
triangular  or  quadrilateral  elements,  or  a  three-dimensional  (3D)  volume  mesh  using  a 
combination  of  hexahedral,  tetrahedral  or  prism  elements. 

Once  the  fluid  domain  has  been  meshed  the  governing  equations  (in  integral  form)  for  the 
conservation  of  mass,  momentum,  energy  and  any  other  relevant  variables  are  applied  to  each 
discrete  control  volume  and  used  to  construct  a  set  of  non-linear  algebraic  equations  for  the 
discrete  dependent  variables.  Fluent  then  solves  the  complete  set  of  coupled  equations  for  all  the 
control  volumes  on  the  mesh  using  either  a  segregated  solver  or  a  coupled  solver.  In  all  the 
calculations  reported  here  only  the  segregated  solver  has  been  used.  In  this  approach  the 
governing  equations  are  solved  sequentially.  Since  these  equations  are  non-linear  they  first  have  to 
be  linearized.  This  can  be  done  either  implicitly  or  explicitly,  although  when  using  the  segregated 
solution  method  within  Fluent  the  code  automatically  linearises  each  discrete  governing  equation 
implicitly  with  respect  to  that  equation's  dependent  variable.  This  produces  a  scalar  system  of 
equations  containing  only  one  equation  per  computational  cell.  A  point  implicit  (Gauss-Siedel) 
linear  equation  solver  is  then  used  in  conjunction  with  an  algebraic  multigrid  (AMG)  method  to 
solve  the  resultant  scalar  system  of  equations  for  the  dependent  variable  in  each  cell.  Since  the 
equations  are  non-linear  several  iterations  of  the  solution  loop  must  be  performed  before  a 
converged  solution  is  obtained. 

The  segregated  solver  thus  solves  for  a  given  fluid  variable  (for  example  U  -  the  x  component  of 
velocity)  by  considering  all  cells  at  a  single  time.  It  then  solves  for  the  next  fluid  variable  (for 
example  V  -  the  y  component  of  velocity)  by  again  considering  all  cells  at  the  same  time.  Each 
iteration  of  the  solution  loop  thus  consists  of  the  following  steps: 

1.  The  values  of  the  fluid  variables  at  the  cell  centres  are  updated  based  on  the  current  solution 
values. 
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2.  The  U,  V  and  W  velocity  components  of  the  momentum  equation  are  each  solved  in  turn 
using  the  current  values  for  the  pressure  and  the  mass  fluxes  through  each  of  the  cell  faces. 

3.  The  pressure  correction  equation  (a  form  of  the  continuity  equation)  is  then  solved  to  obtain 
the  necessary  corrections  to  the  pressure  and  velocity  fields  so  that  the  continuity  equation  is 
satisfied. 

4.  Where  appropriate,  additional  scalar  equations  (such  as  those  describing  transport  of 
turbulence  quantities)  are  solved. 

5.  A  check  for  convergence  is  made. 

6.  The  above  steps  are  repeated  until  the  convergence  criteria  are  met  (all  residuals  less 
than  10'6) 

4.1  Pressure-Velocity  Coupling 

As  described  above,  an  equation  for  each  component  of  the  momentum  equation  and  then  the 
continuity  equation  are  solved  sequentially.  Once  the  three  components  of  velocity  have  been 
calculated  for  each  cell  using  this  sequential  system  the  velocities  may  not  satisfy  the  continuity 
equation  and  so  a  "Poisson-type"  equation  for  a  pressure  correction  is  derived  from  the  continuity 
equation  and  the  linearized  momentum  equations.  This  pressure  correction  equation  is  then  solved 
to  obtain  the  necessary  corrections  to  the  pressure  and  velocity  fields  such  that  continuity  is 
satisfied. 

Although  the  pressure  variable  appears  in  each  of  the  component  momentum  equations  each  of 
these  equations  is  solved  by  treating  the  relevant  component  of  velocity  as  the  unknown  variable 
and  the  pressure  field  in  the  equation  is  taken  to  be  that  from  the  previous  iteration.  In  this 
sequential  procedure  the  continuity  equation  is  used  as  an  equation  for  the  pressure.  However 
pressure  does  not  appear  explicitly  in  the  continuity  equation  for  incompressible  flows  (which  are 
the  only  flows  considered  in  this  report)  and  so  a  procedure  must  be  devised  to  introduce  pressure 
into  this  equation.  Fluent  provides  methods  based  on  the  SIMPLE  (Semi-Implicit  Method  for 
Pressure-Linked  Equations)  family  of  algorithms  to  do  this  [37]. 

The  basic  SIMPLE  algorithm  uses  a  relationship  between  velocity  and  pressure  corrections  to 
enforce  mass  conservation  and  to  obtain  the  pressure  field.  The  SIMPLEC  algorithm  (SIMPLE- 
Consistent)  is  a  variation  of  the  SIMPLE  algorithm  which  uses  a  more  refined  expression  for  the 
variable  flux  through  each  of  the  cell  faces  and  which  can  accelerate  convergence  in  some 
problems  where  the  pressure-velocity  coupling  is  the  main  deterrent  to  obtaining  a  solution  [38]. 
The  PISO  pressure-velocity  coupling  scheme  (Pressure-Implicit  with  Splitting  of  Operators)  is  also 
part  of  the  SIMPLE  family  of  algorithms  and  is  based  on  a  higher  degree  of  approximation  for  the 
relation  between  the  corrections  for  pressure  and  velocity  [39].  The  PISO  algorithm  takes  a  little 
more  CPU  time  per  solver  iteration  but  it  can  dramatically  decrease  the  number  of  iterations 
required  for  convergence,  especially  for  transient  problems.  The  PISO  algorithm  also  allows  Fluent 
to  obtain  solutions  on  highly  skewed  meshes  in  approximately  the  same  number  of  iterations  as 
required  for  more  orthogonal  meshes.  Since  all  the  simulation  results  described  here  were  obtained 
on  fairly  regular  meshes  using  time-independent  calculations  only  the  SIMPLE  algorithm  was 
employed  in  the  calculations. 

4.2  Spatial  Discretization 

Fluent  stores  discrete  values  of  the  variables  at  the  cell  centres,  however  values  of  the  variables  are 
required  at  the  cell  faces  for  the  convection  terms  in  the  equations  and  these  must  be  interpolated 
from  the  cell  centre  values.  This  is  accomplished  using  an  upwind  scheme  and  Fluent  allows  a 
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choice  of  several  different  upwind  schemes:  first-order  upwind,  second-order  upwind,  power  law, 
and  QUICK. 

In  the  first-order  upwind  scheme  quantities  at  cell  faces  are  determined  by  assuming  that  the  cell- 
centre  values  of  any  variable  represent  a  cell-average  value  and  hold  throughout  the  entire  cell. 
Hence  the  face  quantities  are  identical  to  the  cell  quantities.  The  Power-Law  scheme  is  an 
improvement  which  interpolates  the  face  value  of  a  variable  by  using  the  exact  solution  to  a  one¬ 
dimensional  convection-diffusion  equation  [40].  When  the  flow  is  dominated  by  convection  this 
implies  that  the  face  value  of  the  variable  is  effectively  equal  to  the  cell  value  in  the  upwind 
direction.  If  the  flow  is  weak  and  diffusion  stronger  then  the  Power-Law  scheme  amounts  to  a 
simple  linear  average  of  the  value  of  the  variable  at  the  current  cell  location  and  the  upstream  cell. 
The  Second-Order  Upwind  Scheme  provides  true  second  order  accuracy  by  performing  a  Taylor 
series  expansion  of  the  cell-centred  solution  about  the  cell  centroid.  The  QUICK  scheme  also 
allows  calculation  of  a  higher-order  value  of  the  convected  variable  at  the  cell  face  by  using  a 
weighted  average  of  second-order  upwind  and  central  interpolation  [41].  This  scheme  is 
particularly  accurate  on  structured  quadrilateral  and  hexahedral  meshes  where  the  flow  is  aligned 
with  the  grid  direction.  The  Second-Order  Upwind  scheme  was  employed  in  all  of  the  simulations 
described  here. 


5.  Turbulence  Models 


Turbulent  flows  are  characterized  by  velocity  fields  which  fluctuate  rapidly  both  in  space  and 
time.  Since  these  fluctuations  occur  over  several  orders  of  magnitude  it  is  computationally  very 
expensive  to  construct  a  grid  which  directly  simulates  both  the  small  scale  and  high  frequency 
fluctuations  for  problems  of  practical  engineering  significance.  Two  methods  can  be  used  to 
eliminate  the  need  to  resolve  these  small  scales  and  high  frequencies:  Reynolds  Averaging  and 
Filtering. 

In  the  Reynolds  Averaged  approach  all  flow  variables  are  divided  into  a  mean  component  and  a 
rapidly  fluctuating  component  and  then  all  equations  are  time  averaged  to  remove  the  rapidly 
fluctuating  components.  For  the  continuity  equation  the  new  equation  is  identical  to  the  original 
equation,  except  that  the  transported  variables  now  represent  the  mean  flow  quantities.  In  the 
Navier-Stokes  equation  however  new  terms  appear  which  involve  mean  values  of  products  of 
rapidly  varying  quantities.  These  new  terms  are  known  as  the  Reynolds  Stresses,  and  solution  of 
the  Reynolds  Averaged  Navier-Stokes  (RANS)  equation  initially  involves  the  construction  of 
suitable  models  to  represent  these  Reynolds  Stresses.  One  approach  to  this  problem  is  to  treat  the 
time  averaged  terms  as  additional  viscous  stresses  produced  by  the  turbulence  in  the  flow.  In  the 
Boussinesq  approach  the  Reynolds  Stresses  are  assumed  to  have  a  form  identical  to  the  viscous 
stresses  in  the  momentum  equation,  apart  from  a  multiplicative  term  known  as  the  turbulent 
viscosity,  |lit  .  Note  that  this  approach  assumes  that  the  Reynolds  Stresses  are  isotropic,  which  is 
known  to  be  untrue  in  many  cases.  The  problem  then  reduces  to  finding  an  expression  for  pT . 

Fluent  provides  several  turbulence  models  based  on  the  Boussinesq  approach:  the  Spalart- 
Allmaras  model,  the  k-s  model,  and  the  k- co  model.  The  Spalart-Allmaras  model  is  a  relatively 
simple  one-equation  model  which  involves  solving  an  additional  model  transport  equation  for  pT . 
Both  the  k-s  model  and  the  k- co  models  are  inherently  more  complicated  and  involve  finding 
solutions  to  two  additional  model  transport  equations,  one  for  the  turbulent  kinetic  energy  k,  and 
one  for  the  rate  of  dissipation  of  the  tubulent  kinetic  energy  s,  or  in  the  case  of  the  k- co  model,  an 
equation  for  co  (where  co  is  defined  by  co  =  s  /  k).  The  turbulent  viscosity  |lit  is  then  calculated  from 
an  expression  involving  k  and  s  for  the  k-s  model,  or  k  and  co  for  the  k- co  model.  There  are  in  fact 
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three  different  versions  of  the  k-z  model  in  Fluent  and  two  different  versions  of  the  k- co  model,  and 
these  will  be  discussed  in  more  detail  in  the  subsequent  sections. 

A  turbulence  model  which  avoids  making  the  isotropic  Boussinesq  approximation  is  the  Reynolds 
Stress  Model  (RSM).  This  is  the  most  elaborate  turbulence  model  that  Fluent  provides  and  finds  a 
solution  to  the  Reynolds-Averaged  Navier-Stokes  equation  by  solving  additional  transport 
equations  for  each  of  the  individual  Reynolds  stresses,  as  well  as  an  equation  for  the  dissipation 
rate.  This  means  that  four  additional  transport  equations  are  required  in  2D  flows  and  seven 
additional  transport  equations  must  be  solved  in  3D  flows.  Since  the  RSM  accounts  for  the  effects 
of  streamline  curvature,  swirl,  rotation,  and  rapid  changes  in  strain  rate  in  a  more  rigorous  manner 
than  one-equation  (Spalart-Allmaras)  and  two-equation  ( k-z  and  k- 00)  models  it  has  greater 
potential  to  give  accurate  predictions  for  complex  flows.  The  accuracy  of  the  RSM  predictions  is 
still  limited,  however,  by  the  closure  assumptions  used  to  model  various  terms  in  the  exact 
transport  equations  for  the  Reynolds  stresses.  For  this  reason  the  RSM  does  not  always  provide 
results  which  are  superior  to  those  of  simpler  models  for  all  classes  of  flows.  It  is,  of  course, 
necessary  to  use  the  RSM  when  the  flow  features  of  interest  are  the  result  of  anisotropy  in  the 
Reynolds  stresses. 

An  alternative  approach  to  Reynolds  averaging  is  filtering.  The  idea  behind  this  approach  is  to 
filter  the  time-dependent  Navier-Stokes  equation  in  either  Fourier  (wave-number)  space  or 
configuration  (physical)  space.  This  filtering  process  effectively  filters  out  turbulent  eddies  whose 
scales  are  smaller  than  the  filter  width,  which  is  usually  taken  to  be  the  mesh  size.  As  with 
Reynolds  averaging  however,  the  filtering  process  creates  additional  unknown  terms  which  must 
be  modelled  in  order  to  provide  closure  to  the  set  of  equations.  This  approach  is  known  as  Large 
Eddy  Simulation  (LES)  because  the  fluctuations  of  the  large  scale  eddies  (those  having  a  size 
comparable  to  the  main  geometry  of  the  flow)  are  numerically  resolved,  rather  than  being 
averaged  out  as  in  the  RANS  approach.  The  attraction  of  LES  is  that,  by  modelling  less  of  the 
turbulence  (and  solving  more),  the  error  introduced  by  the  turbulence  model  will  be  reduced. 
Fluent  provides  two  methods  to  model  the  subgrid-scale  stresses  resulting  from  the  filtering 
operation:  the  Smagorinsky-Lilly  model  and  the  RNG  (ReNormalization  Group)  subgrid-scale 
model.  All  LES  simulations  require  a  lengthy  time-dependent  run  so  that  statistics  of  the  mean 
flow  quantities  can  be  gathered.  LES  simulations  also  require  a  relatively  fine  grid,  and  so  the 
computational  cost  of  LES  simulations  can  be  quite  prohibitive. 


6.  Computational  Grid 

The  degree  of  resolution  required  of  the  computational  grid  depends  to  some  extent  on  the  choice 
of  turbulence  model  to  be  used  in  the  simulation.  The  k-z  models,  RSM,  and  LES  models  are 
primarily  valid  for  turbulent  core  flows,  i.e.  for  flow  in  the  regions  somewhat  far  from  walls,  while 
the  Spalart-Allmaras  and  k- 00  models  were  designed  to  be  applied  throughout  the  boundary  layer, 
provided  that  the  near-wall  mesh  resolution  is  sufficient.  The  k-z  models  can  still  be  applied  to  wall 
bounded  flows,  however,  by  using  the  concept  of  wall  functions.  In  this  approach  use  is  made  of 
the  universal  behaviour  of  equilibrium  boundary  layer  flows.  It  is  well  known  that  the  near- wall 
region  in  an  equilibrium  boundary  layer  can  be  divided  into  several  distinct  regions  [18].  Very 
close  to  the  wall  the  flow  is  almost  laminar  and  the  molecular  viscosity  plays  a  dominant  role  in 
momentum  transfer.  This  region  is  known  as  the  "viscous  sublayer".  At  much  greater  distances 
from  the  wall,  but  still  well  within  the  boundary  layer,  molecular  viscosity  plays  no  part  and  the 
velocity  profile  is  determined  purely  by  the  turbulent  viscosity.  This  is  known  as  the  outer  layer,  or 
fully-turbulent  layer.  In  between  the  viscous  sublayer  and  the  fully  turbulent  layer  there  is  an 
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interim  layer  where  the  effects  of  molecular  viscosity  and  turbulence  are  equally  important.  This  is 
known  as  the  buffer  layer  or  blending  region. 

When  standard  wall  functions  are  used  the  viscosity  affected  inner  region  (viscous  sublayer  and 
buffer  layer)  is  not  resolved.  Instead,  semi-empirical  formulas  are  used  to  bridge  the  viscosity 
affected  region  between  the  wall  and  the  fully  turbulent  region.  The  use  of  wall  functions  thus 
obviates  the  need  to  modify  the  turbulence  model  to  account  for  the  presence  of  the  wall.  In 
practice,  this  means  that  the  centre  of  the  cell  closest  to  the  wall  must  lie  above  a  certain  height.  If 
y+  denotes  the  (scaled)  co-ordinate  direction  normal  to  a  solid  wall  and  P  denotes  the  centre  point 
of  the  cell  closest  to  the  wall,  then  y+p  should  lie  in  the  range  30<  y+ p<100.  If  this  criterion  is  satisfied 
then  the  boundary  conditions  can  be  satisfied  at  the  point  P  by  using  the  universal  'Tog-law"  for 
the  mean  velocity.  Appropriate  values  for  the  other  variables  can  also  be  derived  at  this  location 
from  the  universal  nature  of  the  flow  in  this  region. 

With  the  above  considerations  in  mind  the  first  computational  mesh  consisted  of  a  structured 
hexahedral  mesh  in  which  the  cell  centres  closest  to  the  wall  were  on  average  2  mm  above  the  flat 
plate  and  the  surface  of  the  wing.  A  total  of  1.6  million  hexahedral  cells  were  used.  This  mesh  was 
used  to  run  simulations  using  the  three  different  k-z  models  available  in  Fluent:  original  k-z,  RNG 
k-z,  and  renormalizeable  k-z,  all  using  the  standard  wall  functions.  Because  the  resolution  of  the 
mesh  was  somewhat  coarser  than  the  resolution  of  the  measurements  used  to  characterize  the 
input  flat  plate  boundary  layer  none  of  these  simulations  used  these  measured  values  as  input 
data.  For  each  of  these  runs  the  y+  values  varied  between  40  and  60  over  the  surface  of  the  plate 
and  between  40  and  100  over  the  leading  edge  of  the  wing.  Each  of  the  simulations  was  able  to 
accurately  calculate  the  pressure  distribution  over  the  surface  of  the  wing,  but  none  of  the 
simulations  showed  the  presence  of  either  a  vortex  or  any  reversed  flow. 

Since  the  first  mesh  proved  to  be  too  coarse  a  finer  mesh  was  constructed  so  that  the  boundary 
layer  profile  measured  upstream  of  the  wing  could  be  used  as  part  of  the  inlet  boundary 
conditions  in  the  simulation.  Because  the  aerofoil  section  had  its  major  axis  aligned  with  the  free- 
stream  flow  a  symmetry  plane  was  created  along  the  chord  of  the  wing  and  only  one  half  of  the 
test  section  was  meshed. 

Devenport  et  al.  [4]  measured  inflow  conditions  1.31  m  upstream  of  the  leading  edge  of  the  wing 
(x/T  =  -18.24).  These  consisted  of  velocity  profiles  measured  at  29  locations  along  the  Y  axis 

between  Y/T  =  0.00354  and  Y/T  =  0.87  at  Z  =  0.0.  The  measured  quantities  were  U/Uref  and  mean 
values  of  u2/  U?ef  /  WV  U?ef  and  uw/  U?ef  •  These  measurements  showed  that  the  boundary  layer 
had  an  initial  thickness  of  approximately  2  cm  and  a  structure  very  similar  to  that  of  a  two- 
dimensional  boundary  layer  under  zero  pressure  gradient  conditions.  The  height  of  the  first  data 
point  above  the  flat  plate  was  0.2  mm  and  the  average  spacing  between  data  points  was  0.8  mm. 

The  computational  grid  extended  from  X  =  -  1.3  m  to  X  =  1.3  m  in  the  X  direction  (X/T  =  -18.13  to 
X/T  =  18.13) ,  Y  =  0.0  to  Y  =  0.25  m  in  the  Y  direction  (Y /T  =  3.49)  and  Z  =  0.0  to  Z  =  0.5  m  in  the  Z 
direction  (Z/T  =  6.97).  The  Y  coordinate  was  divided  into  100  intervals  with  the  first  grid  point 
located  at  a  height  of  approximately  0.1  mm  above  the  flat  plate.  The  size  of  successive  intervals 
increased  with  a  ratio  of  1.04  so  that  the  last  interval  had  a  height  of  approximately  1.0  cm.  The 
surface  of  the  wing  along  the  X  direction  was  meshed  using  110  intervals.  The  first  20  intervals 
were  of  equal  length  (2.1  mm)  and  covered  the  elliptical  part  of  the  nose  while  the  remaining  90 
elements  were  of  varying  length  (average  length  of  3.4  mm)  and  covered  the  NACA  0020  tail  of  the 
wing.  The  first  20  cm  of  the  X  coordinate  upstream  of  the  wing  was  meshed  using  elements  of  2.15 
mm  length  and  similar  griding  was  used  immediately  downstream  of  the  wing.  At  more  distant 
locations  grading  was  employed  so  that  the  length  of  the  elements  increased  slowly  to 
approximately  1  cm  at  the  edges  of  the  mesh.  The  mesh  contained  a  total  of  5  million  hexahedral 
elements.  The  experimentally  measured  values  of  the  velocity  profile  in  the  incoming  boundary 
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layer  flow  (obtainable  in  digital  form  from  [3])  were  used  for  the  velocity  inlet  boundary  condition 
by  constructing  a  Fluent  Profile  file  which  was  read  from  the  Velocity  Inlet  boundary  panel. 


7.  Simulation  Results 


As  well  as  the  standard  wall  functions  described  in  the  previous  section  Fluent  also  provides  both 
non-equilibrium  wall  functions  and  enhanced  wall  functions.  The  two  key  ingredients  of  the  non¬ 
equilibrium  wall  functions  are:  (1)  the  log-law  for  the  mean  velocity  is  sensitised  to  pressure- 
gradient  effects,  (2)  a  two-layer  zonal  model  is  used  to  compute  the  budget  of  turbulent  kinetic 
energy  in  the  wall-neighbouring  cells.  In  a  standard  wall  function  it  is  assumed  that  the  cell- 
averaged  rate  of  production  of  turbulent  kinetic  energy  is  identical  to  the  cell-averaged  dissipation 
rate  of  kinetic  energy.  The  two-layer  zonal  model  employed  in  the  non-equilibrium  wall  functions 
relaxes  this  assumption  by  dividing  each  of  the  wall-neighbouring  cells  into  a  viscous  sublayer  and 
a  fully  turbulent  layer  and  then  computing  the  production  and  dissipation  rates  in  each  of  the 
sublayers.  Non-equilibrium  wall  functions  therefore  extend  the  applicability  of  the  wall  function 
approach  by  including  the  effects  of  pressure  gradients  and  non-equilibrium  kinetic  energy 
production  and  can  give  improvements  in  the  prediction  of  wall  shear  stress  and  heat  transfer  in 
flows  where  such  effects  occur. 

Both  standard  wall  functions  and  non-equilibrium  wall  functions  become  inappropriate  when  the 
flow  contains  either  strong  pressure  gradient  effects  which  result  in  boundary  layer  separation,  or 
strongly  skewed  three-dimensional  boundary  layers.  In  this  situation  a  near-wall  modelling 
approach  combined  with  adequate  mesh  resolution  in  the  near-wall  region  must  be  used.  Fluent 
provides  an  enhanced  wall  treatment  for  these  situations  and  this  approach  can  be  used  with  all 
three  k-z  models  as  well  as  the  RSM  model.  The  enhanced  wall  treatment  is  a  near-wall  modelling 
method  that  combines  a  two-layer  model  with  enhanced  wall  functions.  In  this  approach  the  near¬ 
wall  cells  are  divided  into  a  viscosity-affected  region  and  a  fully-turbulent  region.  The 
demarcation  between  the  two  regions  is  determined  by  a  wall  distance  based  turbulent  Reynolds 
number.  In  the  fully  turbulent  region  either  the  k-z  model  or  the  RSM  model  is  employed,  while  in 
the  viscosity-affected  near  wall  region  the  one-equation  model  of  Wolfstein  [19]  is  used.  In  this 
method  the  momentum  equations  and  the  k  equation  is  retained  but  the  turbulent  viscosity  is 
computed  from  a  simple  expression  involving  k  and  an  appropriate  length  scale.  The  exact  details 
can  be  found  in  the  Fluent  6.0  User's  Guide  Volume  2  [1].  When  using  the  enhanced  wall  treatment 

with  either  the  k-z  models  or  the  RSM  model  it  is  recommended  that  the  y+  value  at  the  wall 
adjacent  cells  should  satisfy  y+  =1,  although  higher  y+  values  are  acceptable  as  long  as  y+  <  5. 

7.1  k-s  Model 

The  k-s  model  is  the  most  widely  used  complete  turbulence  model  and  is  incorporated  in  most 
commercial  CFD  codes.  Fluent  provides  both  the  original  k-z  model  as  developed  by  Jones  and 
Launder  [20]  and  Launder  and  Sharma  [21],  as  well  as  the  RNG  k-z  model  of  Yakhot  and  Orszag 
[22]  which  is  based  on  the  application  of  the  renormalization  group  method  to  the  Navier-Stokes 
equation,  and  the  realizable  k-z  model  of  Shih  et  al.  [23]  which  is  based  on  the  application  of  certain 
constraints  to  the  normal  stresses  to  ensure  that  these  remain  positive  at  all  times. 

The  standard  k-z  model  is  a  semi-empirical  model  based  on  transport  equations  for  the  turbulent 
kinetic  energy  k  and  its  dissipation  rate  s.  The  transport  equation  for  k  was  derived  from  the  exact 
transport  equation  while  the  equation  for  s  was  obtained  using  physical  reasoning  and  cannot  be 
justified  on  an  exact  mathematical  basis.  The  turbulent  viscosity  is  then  specified  by  the  equation 
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vt  =  fc2/ 8,  where  =0.09  is  one  of  five  model  constants.  The  RNG  k-z  model  is  derived  from 
the  instantaneous  Navier-Stokes  equation  using  the  mathematical  technique  known  as  the 
renormalization  group  method.  This  results  in  a  model  with  constants  different  from  those  in  the 
standard  k-z  model  and  additional  terms  in  the  transport  equation  for  8  . 

The  realizable  k-z  model  of  Shih  et  al.  [23]  was  designed  to  remove  two  significant  failings  of 
previous  k-z  models,  these  being  the  tendency  for  the  mean  normal  stresses  (by  definition  these  are 
positive  quantities)  to  become  negative  in  highly  strained  flows,  and  the  rather  empirical  nature  of 
the  model  transport  equation  for  the  dissipation  rate  8.  They  achieved  this  by  introducing  a  new 
eddy-viscosity  formula  in  which  is  no  longer  a  constant  but  has  a  variable  value  which  is  a 
function  of  various  flow  parameters,  and  by  introducing  a  new  model  equation  for  dissipation 
based  on  the  dynamic  equation  for  the  mean-square  vorticity  fluctuation.  Hence  in  the  realizable  k- 
8  model  the  form  of  the  k  equation  is  identical  to  those  in  both  the  standard  k-z  model  and  the  RNG 
k-z  model,  except  for  the  model  constants,  while  the  form  of  the  s  equation  is  quite  different  from 
those  in  the  standard  and  RNG  k-z  models. 

Since  it  is  well  known  that  the  realizable  k-z  model  consistently  outperforms  both  the  standard  k-z 
model  and  the  RNG  k-z  model  in  anything  but  the  simplest  flow  geometries  [18],  only  the 
realizable  k-z  model  has  been  used  for  the  simulation  results  described  here.  This  version  was  run 
on  the  grid  described  in  the  previous  section  using  enhanced  wall  functions  and  the  second  order 
Upwind  scheme  for  the  solution  of  the  momentum,  k  and  8  equations.  The  pressure-velocity 
coupling  used  the  SIMPLE  scheme  and  the  boundaries  simulating  the  wind  tunnel  walls  were 

defined  as  solid  walls  with  specified  zero  shear  stress.  The  y+  values  along  the  centre  line  of  the 
wind  tunnel  floor  upstream  of  the  wing  satisfied  y+  =  3  while  downstream  of  the  wing  they 

satisfied  y+=  4  and  away  from  the  centre  line  approached  y+  =  6.  Along  the  wing  surface  the  y+ 
values  varied  from  approximately  8  on  the  nose  to  4  along  the  trailing  edge.  These  values  are 
reasonably  consistent  with  those  specified  above  and  described  in  the  Fluent  6.0  User's  Guide 

Volume  2  [1].  The  model  constants  used  for  this  particular  implementation  are:  Cis  =  1.44,  C2  =  1.9, 
Gk  =  1.0  and  a8  =  1.2. 

Figures  3  and  4  show  the  pressure  coefficient  around  the  wing  surface,  from  the  leading  edge  to 
the  trailing  edge,  at  two  horizontal  sections,  i.e.  Y/T  =  0.133  and  Y/T  =  1.726.  An  interesting 
feature  of  these  plots  is  that  agreement  with  experiment  is  better  closer  to  the  wall  (Fig.3)  than 
away  from  it  (Fig.  4).  Apsley  and  Leschziner  [10]  found  exactly  the  same  behaviour  in  their 
simulations  and  also  showed  that  this  result  was  not  due  to  grid  dependent  solutions.  They 
concluded  that  the  most  likely  cause  of  the  difference  between  the  simulated  and  experimental 
peak  pressure  shown  in  Fig.  4  was  a  discrepancy  between  the  prescribed  and  experimental  flow 
conditions  at  the  upper  boundary  of  the  computational  domain.  In  the  real  flow  there  was  a  gap 
between  the  top  of  the  wing  and  the  upper  tunnel  wall  while  in  the  computation  (both  in  Apsley 
and  Leschziner's  and  those  reported  here)  the  wing  extended  to  touch  the  upper  surface  of  the 
wind  tunnel.  In  our  view  this  explanation  is  unlikely  as  the  pressure  profile  at  Y/T  =  1.46  is 
identical  to  the  value  at  Y/T  =  1.726,  and  the  value  at  Y/T  =  0.400  also  shows  a  peak  at  this  point. 
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Figure  3.  Pressure  coefficient  along  the  wing  surface  at  a  height  of  Y/T  =  0.133  for  the  Realizable  k-s 
model. 


Figure  4.  Pressure  coefficient  along  the  wing  surface  at  a  height  of  Y/T  =  1.726  for  the  Realizable  k -s 
model. 

Figures  5  through  9  show  a  comparison  between  the  simulated  and  experimental  values  of  U/Uref 
versus  Y/T  at  selected  X/T  values.  For  X/T  values  removed  from  the  centre  of  the  vortex,  i.e. 
those  in  the  range  between  -0.860  and  -0.404,  there  is  fairly  good  agreement  between  the 
simulated  and  experimental  values.  Even  at  the  station  closest  to  the  nose  of  the  wing,  at  X/T  =  - 
0.050,  there  is  still  quite  good  agreement  between  the  experimental  values  and  the  simulated 
values  over  most  of  the  range  of  Y/T.  At  X/T  =  -0.103  there  is  still  reasonable  agreement,  but 
from  X/T  =  -0.154  up  to  X/T  =  -0.351  there  is  significant  disagreement  between  simulation  and 
experiment.  The  Realizable  fc-s  model  essentially  underestimates  the  strength  of  the  reverse  flow  in 
the  vortex  region,  and  this  is  the  same  result  as  noted  by  Apsley  and  Leschziner  [10]  when  using 
the  standard  fc-s  model. 
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Figure  5.  Comparison  of  simulated  and  experimental  values  of  U/Uref  versus  Y/T  at  X/T  =  -0.050  and 
X/T  =  -0.860 for  the  Realizable  k-s  model. 


Figure  6.  Comparison  of  simulated  and  experimental  values  of  U/Uref  versus  Y/T  at  X/T  =  -0.103  and 
X/T  =  -0.667 for  the  Realizable  k-s  model. 


Figure  7.  Comparison  of  simulated  and  experimental  values  of  U/Uref  versus  Y/T  at  X/T  =  -0.154  and 
X/T  =  -0.461  for  the  Realizable  k-s  model. 
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Figure  8.  Comparison  of  simulated  and  experimental  values  of  U/Uref  versus  Y/T  at  X/T  =  -0.251  and 
X/T  =  -0.404 for  the  Realizable  k-s  model. 


Figure  9.  Comparison  of  simulated  and  experimental  values  of  U/Uref  versus  Y/T  at  X/T  =  -0.300  and 
X/T  =  -0.351  for  the  Realizable  k-s  model. 

Figures  10  and  11  show  the  experimental  and  calculated  values  for  V/Uref  versus  Y/T  at  all  the 
X/T  stations.  The  simulated  results  clearly  show  the  same  trend  as  the  experimental  values,  with 
those  X/T  locations  closest  to  the  wing  displaying  the  largest  negative  values  of  velocity. 
Noticeably  absent  however  are  the  positive  components  of  the  velocity  for  those  stations  closest  to 
the  centre  of  the  vortex,  those  at  X/T  =  -0.204  and  X/T  =  -0.251  for  example,  where  the 
experimental  results  show  positive  velocity  components  above  Y/T  =  0.1,  while  the  simulated 
results  are  unable  to  capture  this  behaviour  to  any  extent.  A  more  detailed  comparison  between 
the  experimental  and  simulated  values  is  shown  in  Figures  12  and  13.  These  show  that  at  some  of 
the  X/T  stations  there  is  quite  reasonable  agreement.  For  example,  at  X/T  =  -0.050,  -0.204,  -0.251,  - 
0.300,  -0.404,  and  -0.461  the  simulated  curves  follow  the  experimental  curves  quite  closely,  while 
at  X/T  =  -0.103  and  -0.154  there  is  considerable  disagreement. 

Figure  14  shows  a  comparison  between  the  simulated  and  experimental  values  of  the  turbulent 
kinetic  energy  as  a  function  of  Y/T  for  selected  X/T  values.  Experimentally  the  peak  value  of  the 
turbulent  kinetic  energy  occurs  close  to  the  centre  of  the  vortex  at  X/T  =  -0.351  and  Y/T  =  0.05, 
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with  the  values  either  side  of  the  vortex,  at  X/T  =  -0.050  and  X/T  =  -0.860,  being  considerably 
lower  than  those  near  the  centre.  The  simulated  results  display  quite  a  different  trend.  The  kinetic 
energy  increases  continuously  as  the  distance  to  the  wing  decreases,  and  the  maximum  value 
occurs  not  at  the  vortex  centre  but  rather  at  the  closest  point  to  the  wing  surface.  This  behaviour 
also  occurs  in  the  simulations  of  Apsley  and  Leschziner  [10],  Parneix  et  al.  [11],  and  Rodi  et  al.  [12], 
and  is  a  known  fault  with  all  k-z  models. 


0.1 


-0.35 


Figure  10.  Experimental  values  of  V/Uref  versus  Y/T  at  all  X/T  stations. 


-0.35 


Figure  11.  Simulated  values  of  V/Uref  versus  Y/T  at  all  X/T  stations  for  the  Realizable  k-s  model. 
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—•—V/Uref  X/T  =  -0.050  calc 
V/Uref  X/T  =  -0.050  expt 
—•—V/Uref  X/T  =  -0.103  calc 
-O-  V/Uref  X/T  =  -0.103  expt 
—•—V/Uref  X/T  =  -0.204  calc 
V/Uref  X/T  =  -0.204  expt 
—♦—V/Uref  X/T  =  -0.404  calc 
V/Uref  X/T  =  -0.404  expt 


Figure  12.  Comparison  of  simulated  and  experimental  values  of  V/Uref  versus  Y/T  at  selected  X/T  values  for 
the  Realizable  k-s  model. 


0.07 


—•—V/Uref  X/T  =  -0.251  calc 
—0— V/Uref  X/T  =  -0.251  expt 


-0.18 


Figure  13.  Comparison  of  simulated  and  experimental  values  of  V/Uref  versus  Y/T  at  selected  X/T  values  for 
the  Realizable  k-s  model. 
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Figure  14.  Comparison  of  simulated  and  experimental  values  of  k/U2ref  versus  Y/T  at  selected  X/T  values 
for  the  Realizable  k-s  model. 


7.2  Reynolds  Stress  Model 

This  is  the  most  elaborate  RANS  turbulence  model  provided  in  Fluent  and  since  it  avoids  making 
the  isotropic  Boussinesq  approximation  by  solving  additional  transport  equations  for  each  of  the 
individual  Reynolds  stresses  it  is  expected  to  provide  superior  results  to  any  of  the  other  models. 

The  model  constants  used  for  this  particular  implementation  are:  Ci  =  3.4,  C*i  =  1.8,  C2  =  4.2, 
C3  =  0.8,  C*3  =  1.3,  C4  =  1.25  and  C5  =  0.4. 

Figures  15  through  19  show  a  comparison  between  the  simulated  and  experimental  values  of 
U/Uref  versus  Y/T  at  selected  X/T  values.The  level  of  agreement  between  simulation  and 
experiment  at  X/T  locations  -0.050,  -0.103,  -0.154,  -0.461,  -0.667  and  -0.860  is  very  similar  to  that 
found  for  the  Realizable  k-s  model.  At  locations  close  to  the  vortex  however,  for  X/T  values  of 
-0.251,  -0.300  and  -0.351,  the  RSM  model  provides  considerably  better  agreement  than  the 
Realizable  k-s  model.  The  only  location  at  which  the  Realizable  ks  model  provides  better 
agreement  than  the  RSM  model  occurs  at  X/T  =  -0.404. 

Figures  20  through  22  compare  the  experimental  and  calculated  values  for  V /Uref  versus  Y/T  at  all 
the  X/T  stations.  Comparison  with  Figures  12  and  13  show  that  the  RSM  model  provides 
significant  improvement  over  results  obtained  using  the  Realizable  k-s  model  for  all  of  the  V/Uref 
plots.  There  is  now  much  better  agreement  between  the  simulated  and  experimental  values  at 
X/T  =  -0.103  and  -0.154,  while  at  locations  close  to  the  centre  of  the  vortex,  such  as  X/T  =  -0.251 
and  X/T  =  -0.300  for  example,  the  RSM  model  captures  the  positive  components  of  the  velocity 
which  occur  above  Y/T  =  0.1. 

Figure  23  shows  a  comparison  between  the  simulated  and  experimental  values  of  the  turbulent 
kinetic  energy  as  a  function  of  Y/T  for  selected  X/T  values.  The  results  are  similar  to  those  shown 
in  Figure  14  for  the  Realizable  k-s  model  in  that  the  peak  in  the  kinetic  energy  again  occurs  at  the 
location  closest  to  the  wall,  but  at  least  the  level  of  agreement  between  the  simulated  and 
experimental  value  at  X/T  =  -0.351  is  considerably  improved.  In  summary,  the  RSM  model  seems 
to  provide  better  agreement  with  experiment  for  the  mean  velocity  components  than  does  the 
Realizable  k-s  model,  but  the  turbulent  kinetic  energy  is  still  poorly  predicted. 
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Figure  15.  Comparison  of  simulated  and  experimental  values  of  U/Uref  versus  Y/T  at  X/T  =  -0.050  and 
X/T  =  -0.860 for  the  RSM  model. 


Figure  16.  Comparison  of  simulated  and  experimental  values  of  U/Uref  versus  Y/T  at  X/T  =  -0.103  and 
X/T  =  -0.667 for  the  RSM  model. 


Figure  17.  Comparison  of  simulated  and  experimental  values  of  U/Uref  versus  Y/T  at  X/T  =  -0.154  and 
X/T  =  -0.461  for  the  RSM  model. 
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Figure  18.  Comparison  of  simulated  and  experimental  values  of  U/Uref  versus  Y/T  at  X/T  =  -0.251  and 
X/T  =  -0.404 for  the  RSM  model. 


Figure  19.  Comparison  of  simulated  and  experimental  values  of  U/Uref  versus  Y/T  at  X/T  =  -0.300  and 
X/T  =  -0.351  for  the  RSM  model. 


v/iM  m= 
^v/im  m= 

-•-V/Uref  XTT  = 
9  VIUrefm  = 
-t-V/U-ef  XTT  = 
-^V/U-ef  XfT  = 
— V/U'ef  XfT  = 
-A-V/U-ef  m  = 


-0.050  calc 
-0.050  expt 
-0.103  calc 
-0.103  expt 
-0.204  calc 
-0.204  expt 
■0.404  calc 
■0.404  expt 


Figure  20.  Comparison  of  simulated  and  experimental  values  of  V/Uref  versus  Y/T  at  selected  X/T  values  for 
the  RSM  model 
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Figure  21.  Comparison  of  simulated  and  experimental  values  of  V/Uref  versus  Y/T  at  selected  X/T  values  for 
the  RSM  model 


Figure  22.  Comparison  of  simulated  and  experimental  values  of  V/Uref  versus  Y/T  at  selected  X/T  values  for 
the  RSM  model 


Y/T 


Figure  23.  Comparison  of  simulated  and  experimental  values  of  k/U2ref  versus  Y/T  at  selected  X/T  values 
for  the  RSM. 
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7.3  V2F  Model 

This  model  was  developed  by  Durbin  in  a  number  of  publications  [24-27]  during  the  period  1991- 
1995  and  is  available  in  Fluent  by  payment  of  an  additional  licence  fee.  Conceptually  the  model  lies 
somewhere  between  the  simplicity  of  an  isotropic  k-z  model  and  the  full  anisotropy  provided  by  a 
Reynolds  Stress  model.  It  arose  out  of  a  desire  to  eliminate  the  need  to  use  wall  functions  or  the 
two-layer  approach  to  specify  how  turbulence  varies  near  a  wall.  It  was  developed  from  a 
simplified  Reynolds  Stress  model  and  successfully  includes  the  anisotropy  of  near-wall  turbulence 
while  continuing  to  use  an  eddy  viscosity  assumption.  The  model  is  valid  throughout  the  whole 
flow  domain  and  automatically  becomes  a  near-wall  model  close  to  solid  surfaces.  It  uses  the 

standard  k-z  equations  as  well  as  an  additional  transport  equation  for  the  mean  value  of  v2,  where 
v  represents  the  component  of  the  fluctuating  velocity  normal  to  the  time-averaged  streamline, 
and  an  elliptic  relaxation  equation  for  the  function  f  which  can  be  taken  to  represent  the 
redistribution  of  turbulent  intensity  from  the  streamwise  component  of  the  flow.  The  elliptic 
equation  models  pressure-velocity  correlations  and  accounts  for  the  non-local  effect  of  blocking  by 
the  wall,  ensuring  that  the  velocity  component  normal  to  the  wall  is  inviscidly  brought  to  zero  at 

the  wall.  The  turbulent  viscosity  is  then  given  by  the  expression  vt  =  v  k/ s,  where  has  the 
value  0.2.  The  model  has  a  built-in  realizability  limiter  (Durbin  [28])  and  does  not  require  any 
damping  or  wall  functions  to  correctly  capture  near-wall  turbulence. 

Figures  24  through  28  show  a  comparison  between  the  simulated  and  experimental  values  of 
U/Uref  versus  Y/T  at  selected  X/T  values. The  level  of  agreement  between  simulation  and 
experiment  is  certainly  better  than  that  shown  by  the  Realizable  k-z  model  and  comparable  with 
that  shown  by  the  full  RSM  model.  In  fact  at  X/T  =  -0.050,  -0.103,  -0.154  and  -0.404  there  is  better 
agreement  between  the  simulated  and  experimental  results,  with  the  only  location  showing 
slightly  worse  agreement  being  at  X/T  =  -0.251. 

Figures  29  through  31  compare  the  simulated  and  experimental  values  for  V/Uref  versus  Y/T  at  all 
the  X/T  stations.  Comparison  with  Figures  12  and  13  for  the  Realizable  k-z  model  and  Figures  20 
through  22  for  the  RSM  model  show  that  the  V2F  model  provides  far  better  agreement  with  the 
experimental  values  than  the  Realizable  k-z  model  and  a  similar  degree  of  agreement  as  the  RSM 
model.  At  X/T  =  -0.050  and  -0.103  for  example  the  RSM  model  provides  better  agreement  with 
experiment  than  the  V2F  model,  while  at  X/T  =  -0.204,  -0.251  -0.300  and  -0.351  the  V2F  model 
provides  better  agreement  than  the  RSM  model.  These  figures  show  that  the  V2F  model,  similar  to 
the  RSM  model,  is  capable  of  capturing  the  positive  components  of  the  velocity  which  occur  above 
Y/T  =  0.1  in  the  area  of  the  vortex  centre. 

Figure  32  shows  the  simulated  and  experimental  values  of  the  turbulent  kinetic  energy  as  a 
function  of  Y/T  for  the  same  X/T  locations  shown  in  Figure  14  for  the  Realizable  k-z  model  and 
Figure  23  for  the  RSM  model.  The  results  here  are  an  improvement  on  those  calculated  using  either 
the  Realizable  k-z  model  or  the  RSM  model  close  to  the  wall,  but  agreement  is  worse  near  the 
centre  of  the  vortex.  The  V2F  model  however  still  shows  the  peak  turbulent  kinetic  energy  close  to 
the  wall,  rather  than  at  the  vortex  centre. 
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Figure  24.  Comparison  of  simulated  and  experimental  values  of  U/Uref  versus  Y/T  at  X/T  =  -0.050  and 
X/T  =  -0.860 for  the  V2F  model. 


Figure  25.  Comparison  of  simulated  and  experimental  values  of  U/Uref  versus  Y/T  at  X/T  =  -0.103  and 
X/T  =  -0.667 for  the  V2F  model. 


Figure  26.  Comparison  of  simulated  and  experimental  values  of  U/Uref  versus  Y/T  at  X/T  =  -0.154  and 
X/T  =  -0.461  for  the  V2F  model. 
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Figure  27.  Comparison  of  simulated  and  experimental  values  of  U/Uref  versus  Y/T  at  X/T  =  -0.251  and 
X/T  =  -0.404 for  the  V2F  model. 


Figure  28.  Comparison  of  simulated  and  experimental  values  of  U/Uref  versus  Y/T  at  X/T  =  -0.300  and 
X/T  =  -0.351  for  the  V2F  model. 


Figure  29.  Comparison  of  simulated  and  experimental  values  of  V/Uref  versus  Y/T  at  selected  X/T  values 
for  the  V2F  model. 
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Figure  30.  Comparison  of  simulated  and  experimental  values  of  V/Uref  versus  Y/T  at  selected  X/T  values 
for  the  V2F  model. 


Figure  31.  Comparison  of  simulated  and  experimental  values  of  V/Uref  versus  Y/T  at  selected  X/T  values 
for  the  V2F  model. 
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Figure  32.  Comparison  of  simulated  and  experimental  values  of  k/U2ref  versus  Y/T  at  selected  X/T  values 
for  the  V2F  model. 
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7.4  Spalart- Allmaras  Model 

This  is  the  simplest  turbulence  model  provided  in  Fluent.  It  solves  the  RANS  equations  for  the 
mean  velocity  components  and  pressure  and  uses  the  Boussinesq  hypothesis  to  relate  the  Reynolds 
stresses  to  the  mean  velocity  gradients.  The  turbulent  viscosity  is  then  found  from  the  solution  of 
one  additional  transport  equation.  The  derivation  of  this  equation  is  explained  in  detail  in  the 
original  paper  by  Spalart  and  Allmaras  [29].  Rather  than  being  derived  in  a  more  or  less 
fundamental  way  from  moments  of  the  Navier-Stokes  equation,  as  are  the  equations  for  k  and  s  in 
the  k-z  model,  the  transport  equation  for  the  turbulent  viscosity  is  simply  postulated,  based  on  a 
combination  of  empiricism,  dimensional  analysis,  Galilean  invariance,  and  selective  dependence 
on  the  molecular  viscosity.  In  its  original  formulation  it  contains  production,  destruction  and 
diffusion  terms,  with  both  the  production  and  destruction  terms  being  a  function  of  the  distance  to 
the  wall,  and  the  production  term  also  being  a  function  of  the  magnitude  of  the  vorticity.  The 
model  was  designed  specifically  for  aerospace  applications  involving  wall-bounded  flows  and  has 
been  shown  to  give  good  results  for  boundary  layers  subjected  to  adverse  pressure  gradients. 

Fluent  contains  both  the  original  version  of  the  model  as  well  as  a  modified  form  which  takes  into 
account  the  effect  of  mean  strain  on  the  turbulence  production.  This  modification  combines 
measures  of  both  rotation  and  strain  tensors  in  the  definition  of  the  production  term.  Including 
both  the  rotation  and  strain  tensors  reduces  the  production  of  eddy  viscosity  and  consequently 
reduces  the  eddy  viscosity  itself  in  regions  where  the  measure  of  the  vorticity  exceeds  that  of  the 
strain  rate.  This  is  typically  found  in  vortical  flows  near  the  core  of  a  vortex  subjected  to  pure 
rotation  where  turbulence  is  known  to  be  supressed.  This  is  the  opposite  of  the  case  under  study 
here,  where  the  turbulence  production  in  the  core  of  the  vortex  is  considerably  enhanced,  and  so 
only  the  original  version  of  the  model  has  been  used  in  this  report. 

The  original  Spalart- Allmaras  model  was  effectively  a  low  Reynolds  number  model  and  required 
fine  resolution  in  the  viscosity  affected  region  of  the  boundary  layer.  In  Fluent,  however,  the 
model  has  been  implemented  to  use  wall  functions  when  the  mesh  resolution  is  not  sufficiently 
fine.  Guidelines  for  the  use  of  this  model,  as  described  in  the  Fluent  6.1  Tutorial  Guide  Number  3 

[30],  indicate  that  the  y+  value  of  the  wall  adjacent  cells  should  be  either  very  small  (on  the  order  of 

y+  =  1),  or  approximately  30  or  greater.  Paciorri  et  al.  [31]  used  the  Spalart-Allmaras  model  to 
simulate  flow  around  a  hollow-cylinder  flare  problem  and  found  excellent  agreement  with  the 
experimental  data  provided  that  the  wall  y+  value  did  not  exceed  a  value  of  2.  On  a  more  coarse 
mesh,  where  the  y+  value  exceeded  6,  the  agreement  with  experiment  was  severely  degraded. 

Since  the  y+  values  on  the  present  grid  were  of  the  order  of  5  and  6,  a  new  mesh  was  constructed 
with  finer  resolution  in  the  areas  close  to  the  wall.  The  height  of  the  first  layer  of  cells  above  the 
solid  surfaces  was  changed  from  0.1mm  to  approximately  0.03  mm  and  stronger  growth  ratios 

were  used  so  that  the  cell  sizes  at  the  edges  of  the  mesh  were  approximately  2  cm.  The  resulting  y+ 
values  around  the  base  of  the  wing  varied  between  0.5  to  1.5,  while  over  the  leading  edge  of  the 
wing  the  maximum  y+  value  was  approximately  2.  The  model  constants  used  for  this  particular 
implementation  are:  Cbi  =  0.1335,  Cb2  =  0.622,  av  =  2/3,  Cvi  =  7.1,  Cwi  =  Cbi/iQ  +  (1+Cb2)/  / 

Cw2  =  0.3,  Cw3  =  2.0  and  k  =  0.4187. 

Figures  33  through  37  show  a  comparison  between  the  simulated  and  experimental  values  of 
U/Uref  versus  Y/T  at  selected  X/T  values.  The  level  of  agreement  between  simulation  and 
experiment  is  again  noticeably  better  than  that  shown  by  the  Realizable  k-s  model  and  comparable 
with  that  shown  by  the  full  RSM  model  or  the  V2F  model.  At  X/T  =  -0.050,  -0.103,  -0.154,  -0.461,  - 
0.667  and  -0.860  all  the  models  show  very  similar  levels  of  agreement,  although  at  X/T  =  -0.251,  - 
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0.404  and  -0.461  both  the  V2F  model  and  the  Spalart-Allmaras  model  show  better  agreement  with 
experiment,  and  the  Spalart-Allmaras  results  are  remarkably  similar  to  those  of  the  V2F  model.  At 
X/T  =  -0.300  and  -0.351  however  the  Spalart-Allmaras  results  are  inferior  to  those  of  either  the 
RSM  or  V2F  models,  although  they  are  still  better  than  those  calculated  by  the  Realizable  k-z 
model. 

Figures  38  through  40  compare  the  simulated  and  experimental  values  for  V/Uref  versus  Y/T  at  all 
the  X/T  stations.  These  demonstrate  again  that  the  level  of  agreement  between  simulation  and 
experiment  is  noticeably  better  than  that  shown  by  the  Realizable  k-z  model  and  comparable  with 
that  shown  by  the  full  RSM  model  or  the  V2F  model.  Only  at  X/T  =  -0.251  does  the  Spalart- 
Allmaras  model  fail  to  predict  the  positive  velocity  component  which  occurs  close  to  the  wall  in 
the  area  of  the  vortex  centre.  Note  that  since  the  Spalart-Allmaras  model  solves  a  transport 
equation  for  the  turbulent  viscosity  rather  than  the  turbulent  kinetic  energy  we  are  unable  to  make 
a  comparison  between  simulated  and  experimental  values  of  the  turbulent  kinetic  energy  as  a 
function  of  Y/T  at  the  X/T  locations. 


Figure  33.  Comparison  of  simulated  and  experimental  values  of  U/Uref  versus  Y/T  at  X/T  =  -0.050  and 
X/T  =  -0.860 for  the  Spalart-Allmaras  model. 


Figure  34.  Comparison  of  simulated  and  experimental  values  of  U/Uref  versus  Y/T  at  X/T  =  -0.103  and 
X/T  =  -0.667 for  the  Spalart-Allmaras  model. 
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Figure  35.  Comparison  of  simulated  and  experimental  values  of  U/Uref  versus  Y/T  at  X/T  =  -0.154  and 
X/T  =  -0.461  for  the  Spalart-Allmaras  model. 


Figure  36.  Comparison  of  simulated  and  experimental  values  of  U/Uref  versus  Y/T  at  X/T  =  -0.251  and 
X/T  =  -0.404 for  the  Spalart-Allmaras  model. 


Figure  37.  Comparison  of  simulated  and  experimental  values  of  U/Uref  versus  Y/T  at  X/T  =  -0.300  and 
X/T  =  -0.351  for  the  Spalart-Allmaras  model. 
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Figure  38.  Comparison  of  simulated  and  experimental  values  of  V/Uref  versus  Y/T  at  selected  X/T  values 
for  the  Spalart-Allmaras  model. 


Figure  39.  Comparison  of  simulated  and  experimental  values  of  V/Uref  versus  Y/T  at  selected  X/T  values 
for  the  Spalart-Allmaras  model. 


Figure  40.  Comparison  of  simulated  and  experimental  values  of  V/Uref  versus  Y/T  at  selected  X/T  values 
for  the  Spalart-Allmaras  model. 
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7.5  k-co  Model 

The  k-oo  model  is  the  second  most  widely  used  two  equation  model.  It  was  developed  by  Wilcox  in 
1988  [32],  although  an  improved  version  appeared  in  the  book  by  Wilcox  in  1998  [33].  In  this 
model  the  expression  for  the  turbulent  viscosity  and  the  equation  for  the  turbulent  kinetic  energy 
are  the  same  as  those  used  in  the  k-z  model.  The  difference  between  the  two  models  lies  in  the 
second  equation,  which  Wilcox  writes  as  an  equation  for  the  variable  o,  which  is  the  dissipation 
per  unit  turbulence  kinetic  energy  and  is  defined  simply  as  c o  =  z  /  k.  Wilcox  has  shown  that,  for 
boundary  layer  flows,  the  k-co  model  is  superior  to  the  k-z  model  both  in  its  treatment  of  the 
viscous  near  wall  region  and  in  its  accounting  for  the  effects  of  streamwise  pressure  gradients.  The 
original  model  had  problems  near  non-turbulent  free  stream  boundaries  because  a  non-zero  (non¬ 
physical)  boundary  condition  is  required  for  oo  and  the  calculated  flow  was  found  to  be  sensitive  to 
the  value  specified.  The  modified  k-co  model  described  by  Wilcox  in  1998  [33]  has  removed  this 
problem. 

A  different  version  of  the  k- co  model,  the  shear-stress  transport  (SST)  k-co  model,  was  described  by 
Menter  in  1994  [34].  In  this  model  the  turbulent  viscosity  is  modified  to  account  for  the  transport  of 
the  principal  turbulent  shear  stress.  This  feature  gives  the  model  an  advantage  in  terms  of 
performance  over  both  the  standard  k-co  and  k-z  models.  Other  modifications  include  the  addition 
of  a  cross-diffusion  term  in  the  oo  equation  and  the  use  of  a  blending  function.  Close  to  walls  the 
blending  function  is  zero,  whereas  remote  from  walls  the  blending  function  is  unity  and  the  model 
reverts  to  the  standard  k-z  model.  Hence  the  SST  k-co  model  blends  the  k-co  model,  which  is 
superior  in  the  near-wall  region,  with  the  k-z  model,  which  is  less  sensitive  to  the  level  of 
turbulence  in  the  free  stream. 

Fluent  provides  both  the  original  k-co  model  as  well  as  the  SST  k-co  model  and  also  provides  each 
of  these  models  as  low-Reynolds  number  models  (fine  grids)  as  well  as  high-Reynolds  number 
models  (coarse  grids).  Only  the  low-Reynolds  number  SST  k-co  model  has  been  used  in  this  report 
as  this  represents  the  most  refined  version  of  the  model  to  date.  The  simulation  was  performed 
using  the  grid  described  in  the  previous  section  to  ensure  maximum  resolution  in  the  near  wall 

region.  The  model  constants  used  for  this  particular  implementation  are:  a*oo  =1,  otoo  =0.52,  ao  = 
1/9,  (3*oo  =0.09,  Pf  =0.072,  Rp  =  8,  Rk  =  6,  Rw  =  2.95,  ^  =  1.5,  Mt0  =  0.25,  ak  =  2.0,  =  2.0. 

Figures  41  through  45  show  a  comparison  between  the  simulated  and  experimental  values  of 
U/Uref  versus  Y/T  at  selected  X/T  values.  The  level  of  agreement  between  simulation  and 
experiment  is  again  noticeably  better  than  that  shown  by  the  Realizable  k-z  model  but  perhaps  not 
quite  as  good  as  that  shown  by  the  full  RSM  model.  The  calculated  results  appear  on  average  to 
have  the  same  level  of  agreement  with  experiment  as  those  calculated  by  either  the  V2F  model  or 
the  Spalart-Allmaras  model. 

Figures  46  through  48  compare  the  simulated  and  experimental  values  for  V/Uref  versus  Y/T  and 
the  results  show  a  similar  trend  to  those  shown  by  the  U/Uref  calculations;  the  simulated  results 
are  again  better  than  those  calculated  by  using  the  Realizable  k-z  model  but  not  as  accurate  as  those 
calculated  using  the  full  RSM  model.  The  level  of  accuracy  again  seems  to  be  similar  to  that 
obtained  when  using  either  the  V2F  model  or  the  Spalart-Allmaras  model. 
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Figure  41.  Comparison  of  simulated  and  experimental  values  of  U/Uref  versus  Y/T  at  X/T  =  -0.050  and 
XfT  -  -0.860 for  the  SST  k-co  model 


Figure  42.  Comparison  of  simulated  and  experimental  values  of  U/Uref  versus  Y/T  at  X/T  =  -0.103  and 
X/T  =  -0.667 for  the  SST  k-co  model. 


Figure  43.  Comparison  of  simulated  and  experimental  values  of  U/Uref  versus  Y/T  at  X/T  =  -0.154  and 
X/T  =  -0.461  for  the  SST  k-co  model. 
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Figure  44.  Comparison  of  simulated  and  experimental  values  of  U/Uref  versus  Y/T  at  X/T  =  -0.251  and 
XfT  -  -0.404 for  the  SST  k-co  model. 


Figure  45.  Comparison  of  simulated  and  experimental  values  of  U/Uref  versus  Y/T  at  X/T  =  -0.300  and 
X/T  =  -0.351  for  the  SST  k-co  model. 

Figure  49  shows  the  simulated  and  experimental  values  of  the  turbulent  kinetic  energy  as  a 
function  of  Y/T  for  the  same  X/T  locations  shown  in  Figure  14  for  the  Realizable  k-s  model.  Figure 
23  for  the  RSM  model  and  Figure  32  for  the  V2F  model.  Given  that  the  SST  k-co  model  was 
designed  to  be  integrated  through  the  viscous  sublayer  the  results  here  are  disappointing;  the 
simulated  result  at  X/T  =  -0.050  is  worse  than  that  calculated  by  any  of  the  other  models  and  the 
results  at  X/T  =  -0.351  and  -0.860  are  comparable  to  those  calculated  using  any  of  the  other 
models.  The  simulated  results  still  show  the  same  incorrect  trend,  with  the  peak  turbulent  kinetic 
energy  occuring  closest  to  the  wall,  rather  than  at  the  vortex  centre. 
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Figure  46.  Comparison  of  simulated  and  experimental  values  of  V/Uref  versus  Y/T  at  selected  X/T  values 
for  the  SST  k-co  model  model. 


Figure  47.  Comparison  of  simulated  and  experimental  values  of  V/Uref  versus  Y/T  at  selected  X/T  values 
for  the  SST  k-co  model  model. 
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Figure  48.  Comparison  of  simulated  and  experimental  values  of  V/Uref  versus  Y/T  at  selected  X/T  values 
for  the  SST  k-co  model. 
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Figure  49. 


Comparison  of  simulated  and  experimental  values  of  k/U2ref  versus  Y/T  at  selected  X/T  values 
for  the  SST  k-co  model. 


8.  Grid  Independence 

To  ensure  that  the  results  presented  in  the  previous  sections  were  grid  independent,  simulations 
for  each  of  the  turbulence  models  were  run  on  either  two  or  three  different  grids  having  varying 
degrees  of  refinement.  In  this  section  we  discuss  the  accuracy  of  the  simulation  results  for  each  of 
the  turbulence  models  separately. 

Simulations  using  the  realizable  k-s  model  were  performed  on  both  the  second  mesh  described  in 
Section  6.0  (henceforth  referred  to  as  the  coarse  grid)  and  the  mesh  described  in  Section  7.4,  which 
was  constructed  for  the  Spalart-Allmaras  model  (which  will  be  referred  to  as  the  fine  mesh). 

Figures  50  and  51  show  plots  of  U/Uref  and  V/Uref  as  a  function  of  Y/T  at  several  representative 
X/T  locations  for  simulations  performed  on  both  the  coarse  and  fine  meshes.  Figure  50  shows  that 
the  U/Uref  velocity  profiles  calculated  on  the  two  meshes  are  identical,  even  in  the  highly 
turbulent  region  near  the  vortex  core,  i.e.  at  X/T  =  -0.300.  The  V/Uref  velocity  profiles  shown  in 
Figure  51  however  show  that  grid  independence  has  not  been  completely  achieved.  The  profiles  at 
X/T  =  -0.050  and  X/T  =  -0.103  for  example,  while  very  similar,  show  differences  of  up  to  10%  at 
certain  Y/T  locations.  Strangely  enough  the  profiles  closer  to  the  vortex  core,  for  example  at  X/T  = 
-0.154,  are  again  almost  identical.  To  achieve  complete  grid  independence  for  this  model  would 
probably  require  a  finer  mesh  than  that  described  in  Section  7.4.  Such  a  mesh  would  be  expected  to 
improve  the  accuracy  of  the  simulations  by  only  a  few  percent  however,  and  the  experimental 
values  shown  in  Figure  51  for  X/T  =  -0.050  and  X/T  =  -0.103  show  that  improvements  to  the 
simulated  results  of  this  order  would  have  little  effect  on  the  overall  agreement  between  the 
simulated  and  experimental  results.  We  conclude  that  the  results  presented  in  Section  7.1  for  the 
realizable  k-s  model  show  sufficient  grid  independence  as  to  make  the  results  representative  of  the 
degree  of  agreement  between  simulation  and  experiment  which  can  be  achieved  using  this  model. 
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Figure  50.  Plots  of  U/Uref  versus  Y/T  at  selected  X/T  locations  for  simulations  on  the  coarse  and  fine  meshes 
for  the  realizable  k-s  model. 
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Figure  51.  Plots  of  V/Uref  versus  Y/T  at  selected  X/T  locations  for  simulations  on  the  coarse  and  fine  meshes 
for  the  realizable  k-s  model. 

For  the  V2F  model  simulations  were  performed  on  both  the  coarse  and  fine  grids  and  comparisons 
between  the  results  obtained  on  both  grids  are  shown  in  Figures  52  and  53,  which  show  plots  of 

U/Uref  and  V/Uref  respectively  versus  Y/T  at  selected  X/T  locations.  Figure  52  shows  that  the 

U/Uref  velocity  profiles  calculated  on  the  two  meshes  are  virtually  identical.  There  is  a  very  slight 
difference  in  the  region  of  the  vortex  centre,  near  X/T  =  -0.300,  but  this  only  on  the  order  of  a  few 
percent.  The  V/Uref  velocity  profiles  in  Figure  53  show  a  similar  degree  of  convergence  as  that 
achieved  using  the  realizable  fc-s  model.  Away  from  the  vortex  centre  the  simulated  velocity 
profiles  on  the  two  grids  are  almost  identical,  but  at  X/T  =  -0.300  there  is  again  a  difference  of  a 
few  percent  over  part  of  the  Y/T  range.  Although  the  simulated  profiles  on  the  coarse  and  fine 
meshes  are  not  identical  everywhere  the  differences  are  no  more  than  a  few  percent  in  isolated 
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regions  and  are  certainly  significantly  less  than  the  scatter  in  the  experimental  data  in  these 
regions,  and  so  we  conclude  that  the  results  shown  in  Section  7.3  for  the  V2F  model  are  again 
representative  of  the  degree  of  agreement  between  simulation  and  experiment  which  can  be 
achieved  using  this  model. 

U/Uref  at  X/T  =  -0.050  fine 


Figure  52.  Plots  of  U/Uref  versus  Y/T  at  selected  X/T  locations  for  simulations  on  the  coarse  and  fine  meshes 
for  the  V2F  model. 


Figure  53.  Plots  of  V/UTef  versus  Y/T  at  selected  X/T  locations  for  simulations  on  the  coarse  and  fine  meshes 
for  the  V2F  model. 

Figures  54  and  55  show  plots  of  U/Uref  and  V/Uref  versus  Y/T  at  selected  X/T  locations  for  the 
Spalart-Allmaras  model.  For  this  model  simulations  were  run  on  three  different  grids;  the  coarse 
and  fine  grids  referred  to  previously,  as  well  as  an  even  finer  grid  which  was  obtained  from  the 
fine  grid  using  the  Fluent  grid  adaption  feature.  Using  the  hanging  node  method  cells  marked  for 
refinement  are  refined  by  isotropically  subdividing  each  cell,  so  that  a  hexahedral  cell  for  example 
is  split  into  eight  equivalent  hexahedra.  Cells  can  be  marked  for  refinement  using  a  number  of 
criteria,  but  in  this  case,  knowing  the  location  of  the  vortex,  the  cells  were  refined  in  a  region  which 
completely  encapsulated  the  vortical  flow.  This  effectively  doubled  the  number  of  cells  in  the 
mesh,  from  approximately  1  million  hexahedral  cells  to  approximately  2  million  cells. 
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From  the  plots  shown  in  Figure  54  it  could  be  concluded  that  grid  convergence  has  been  achieved 
on  the  fine  mesh,  but  inspection  of  the  velocity  profiles  in  Figure  55  show  that  there  are  very  small 
differences  between  the  profiles  calculated  on  the  fine  and  adapted  meshes  over  part  of  the  Y/T 
ranges  for  some  of  the  X/T  locations.  These  again  occur  very  close  to  the  vortex  core,  at  X/T  =  - 
0.204  for  example.  Elsewhere  there  is  excellent  agreement  between  the  profiles  calculated  on  the 
three  different  grids.  Considering  the  level  of  agreement  between  the  experimental  results  and  the 

V/Uref  velocity  profiles  shown  in  Figures  38  to  40  however  it  was  considered  that  further  grid 
refinement  in  the  region  of  the  vortex  core  was  unnecessary  as  it  was  unlikely  to  significantly 
increase  the  overall  level  of  agreement  between  the  simulated  and  the  experimental  results  for  the 
Spalart-Allmaras  model. 


-•-U/Uref  at  X/T  = 
o  U/Uref  at  X/T  = 
-•-U/Uref  at  X/T  = 
a  U/Uref  at  X/T  = 
-•-U/Uref  at  X/T  = 
O  U/Uref  at  X/T  = 
o  U/Uref  at  X/T  = 
U/Uref  at  X/T  = 
o  U/Uref  at  X/T  = 


-0.050  coarse 
-0.050  fine 
-0.860  coarse 
-0.860  fine 
-0.404  coarse 
-0.404  fine 
-0.404  adapted 
-0.154  adapted 
-0.154  fine 


-0.6 


Figure  54.  Plots  of  U/Uref  versus  Y/T  at  selected  X/T  locations  for  simulations  on  the  coarse ,  fine  and 
adapted  meshes  for  the  Spalart-Allmaras  model. 


Figure  55.  Plots  of  V/Uref  versus  Y/T  at  selected  X/T  locations  for  simulations  on  the  coarse ,  fine  and 
adapted  meshes  for  the  Spalart-Allmaras  model. 
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For  the  k- co  model  simulations  were  performed  on  the  fine  mesh  and  the  adapted  mesh  described 
in  the  preceeding  paragraphs.  Figure  56  shows  plots  of  U/Uref  versus  Y/T  at  X/T  =  -0.050,  X/T  = 
-0.404  and  X/T  =  -0.860  and  the  profiles  calculated  on  the  two  grids  are  identical.  Figure  57  shows 
plots  of  U/Uref  versus  Y/T  at  X/T  =  -0.154  and  X/T  =  -0.251.  At  X/T  =  -0.154  the  profiles  on  the 
two  different  meshes  are  again  identical,  but  at  X/T  =  -0.251  there  are  differences  of  a  few  percent 
at  the  smallest  Y/T  values.  Examination  of  the  experimental  data  at  this  X/T  station  however 
shows  that  there  is  little  point  in  further  refining  the  grid  in  this  region  as  the  velocity  profile 
simulated  on  this  even  finer  grid  would  be  unlikely  to  provide  significant  improvement  in  the 
agreement  between  the  simulated  and  experimental  data. 


Figure  56.  Plots  of  U/Uref  versus  Y/T  at  selected  X/T  locations  for  simulations  on  the  fine  and  adapted 
meshes  for  the  k-co  model 

Figure  58  for  the  plots  of  V/Uref  versus  Y/T  at  X/T  =  -0.050,  X/T  =  -0.154,  X/T  =  -0.300  and 
X/T  =  -0.404  shows  a  similar  trend.  There  is  excellent  agreement  between  the  velocity  profiles 
calculated  on  the  two  grids  for  X/T  =  -0.050,  X/T  =  -0.154  and  X/T  =  -0.404,  but  at  X/T  =  -0.300 
there  is  a  significant  difference  between  the  profiles  over  a  considerable  part  of  the  Y/T  range. 
Again,  examination  of  the  experimental  data  for  this  X/T  location  shows  that  further  grid 
refinement  in  this  region  is  unlikely  to  significantly  improve  the  agreement  between  the  simulated 
and  experimental  data. 


Figure  57.  Plots  of  U/Uref  versus  Y/T  at  selected  X/T  locations  for  simulations  on  the  fine  and  adapted 
meshes  for  the  k-co  model 
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Figure  58.  Plots  of  V/Uref  versus  Y/T  at  selected  X/T  locations  for  simulations  on  the  fine  and  adapted 
meshes  for  the  k-co  model 

The  simulations  on  the  different  meshes  presented  so  far  in  this  section  have  demonstrated 
effective  grid  convergence  for  both  the  realizable  fc-e  model  and  the  V2F  model  on  the  coarse  mesh, 
and  for  the  Spalart-Allmaras  model  and  the  k-co  model  on  the  fine  mesh.  For  the  RSM  model 
however  the  situation  is  more  complex.  Figures  59  through  66  show  simulations  for  the  RSM 
model  on  both  the  coarse  and  fine  meshes  using  both  non-equilibrium  wall  functions  as  well  as 
advanced  wall  functions.  Both  wall  functions  were  trialed  because  the  Fluent  6.0  User's  Guide 
states  "Using  the  RSM  creates  a  high  degree  of  coupling  between  the  momentum  equations  and 
the  turbulent  stresses  in  the  flow,  and  thus  the  calculation  can  be  more  prone  to  stability  and 
convergence  difficulties  than  with  the  k-e  models."  [35].  The  following  recommendation  is  also 
noted  "When  calculating  with  RSM  it  is  recommended  that  the  Non-Equilibrium  Wall  Functions 
are  used"  [36].  Whilst  the  reason  for  this  recommendation  was  not  given,  it  was  certainly  noted 
that  convergence  of  the  RSM  on  both  the  coarse  and  fine  mesh  was  easily  obtained  using  non¬ 
equilibrium  wall  functions,  while  convergence  on  any  mesh  using  enhanced  wall  functions  was 
extremely  difficult  to  obtain.  Also,  it  should  be  noted  that  much  better  convergence  from  the  point 
of  view  of  residuals  was  obtained  using  non-equilibrium  wall  functions.  On  either  the  coarse  or 
fine  grid  the  residuals  were  observed  to  reduce  by  up  to  six  orders  of  magnitude  using  the  non¬ 
equilibrium  wall  functions,  whilst  when  using  the  enhanced  wall  functions  the  residuals  typically 
reduced  by  only  three  orders  of  magnitude  and  then  the  calculation  became  unstable.  Hence, 
whilst  a  calculation  using  enhanced  wall  functions  on  the  fine  grid  might  be  expected  to  give 
superior  results  to  a  calculation  on  the  coarse  grid  using  non-equilibrium  wall  functions,  this  is  not 
necessarily  the  case  because  of  the  poor  convergence  (from  the  residuals  view  point)  of  the 
calculation  on  the  fine  mesh  using  enhanced  wall  functions. 

Examination  of  Figures  59  through  66  shows  some  interesting  features.  For  many  of  the  X/T 
locations  away  from  the  centre  of  the  vortex,  for  example  at  X/T  =  -0.050,  -0.103,  -0.154  and  -0.860 
there  is  fairly  good  agreement  between  simulations  performed  on  the  coarse  and  fine  grids  using 
the  non-equilibrium  wall  functions  and  those  on  the  fine  grid  using  the  enhanced  wall  functions. 
Even  at  X/T  =  -0.204  and  -0.404  the  agreement  between  these  three  simulations  is  still  reasonable. 

Only  at  X/T  =-0.300  for  U/Uref  and  at  X/T  =  -0.251  for  V/Uref  are  there  significant  differences 
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between  the  three  simulations.  In  both  these  cases  however  it  is  interesting  to  note  that  the 
experimental  results  agree  more  closely  with  the  simulation  performed  on  the  coarse  mesh  using 
the  non-equilibrium  wall  functions.  This  was  also  found  to  be  the  case  at  other  X/T  locations. 
These  simulations  are  obviously  not  fully  converged,  but  pursuing  convergence  on  more  and  more 
refined  grids  using  non-equilibrium  wall  functions  is  questionable  in  any  case,  and  it  is  also  not 
apparent  that  convergence  can  be  obtained  using  enhanced  wall  functions.  For  these  reasons  the 
results  presented  in  Section  7.2  use  the  simulations  performed  on  the  coarse  mesh  using  the  non¬ 
equilibrium  wall  functions. 


Figure  59.  Plots  of  U/Uref  versus  Y/T  at  X/T  =  -0.050  and  X/T  =  -0.860  for  simulations  on  the  coarse  and 
fine  meshes  for  the  RSM  model. 


Figure  60.  Plots  of  U/Uref  versus  Y/T  at  X/T  =  -0.204  for  simulations  on  the  coarse  and  fine  meshes  for  the 
RSM  model. 
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Figure  61.  Plots  of  U/Uref  versus  Y/T  at  X/T  =  -0.300  for  simulations  on  the  coarse  and  fine  meshes  for  the 
RSM  model. 


Figure  62.  Plots  of  U/Uref  versus  Y/T  at  X/T  =  -0.404  for  simulations  on  the  coarse  and  fine  meshes  for  the 
RSM  model. 


Figure  63.  Plots  of  V/Uref  versus  Y/T  at  X/T  =  -0.050  for  simulations  on  the  coarse  and  fine  meshes  for  the 
RSM  model. 
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Figure  64.  Plots  of  V/Uref  versus  Y/T  at  X/T  =  -0.103  for  simulations  on  the  coarse  and  fine  meshes  for  the 
RSM  model. 


Figure  65.  Plots  of  V/Uref  versus  Y/T  at  X/T  =  -0.154  for  simulations  on  the  coarse  and  fine  meshes  for  the 
RSM  model. 


Figure  66.  Plots  of  V/Uref  versus  Y/T  at  X/T  =  -0.251  for  simulations  on  the  coarse  and  fine  meshes  for  the 
RSM  model. 
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9.  Discussion  and  Conclusion 


The  aim  of  this  report  was  to  investigate  and  document  the  accuracy  with  which  the  Fluent  CFD 
code  and  the  various  turbulence  models  contained  within  it  could  simulate  the  horseshoe  vortex 
formed  in  the  Devenport-Simpson  experiment  [4].  Simulations  were  conducted  using  each  of  the 
turbulence  models  and  the  calculated  results  were  compared  with  experimental  results  obtained 
from  an  extensive  database  of  results  available  at  the  University  of  Surrey's  website  [3].  In  all  cases 
comparisons  were  made  in  the  plane  of  symmetry  upstream  of  the  wing  for  the  mean  components 
of  velocity  in  the  direction  of  the  free-stream  and  normal  to  the  flat  plate  at  eleven  stations 
upstream  of  the  wing,  as  well  as  the  pressure  distribution  along  the  surface  of  the  wing  at  selected 
locations. 

All  of  the  calculations  accurately  reproduced  the  pressure  distributions  around  the  wing  and  only 
results  for  the  k-z  model  were  presented  as  the  remaining  models  displayed  identical  results, 
which  was  as  expected.  The  mean  velocity  components  in  the  free-stream  direction  and  normal  to 
the  flat  plate  were  reasonably  well  simulated  by  all  the  turbulence  models,  except  in  the  vicinity  of 
the  vortex  core.  The  realizable  k-z  model  and  the  Spalart-Allmaras  model  were  less  accurate  than 
the  other  models  in  simulating  velocity  components  in  the  vortex  core,  while  the  remaining 
models,  the  RSM,  the  V2F  model  and  the  k-o  model  all  displayed  similar  levels  of  accuracy  with 
regard  to  simulations  of  the  mean  velocity  components  in  the  vortex  core.  None  of  the  models  was 
able  to  accurately  simulate  the  correct  trend  for  the  behaviour  of  the  mean  kinetic  energy  as  a 
function  of  Y/T  at  selected  X/T  values.  Experimentally  the  peak  value  of  the  turbulent  kinetic 
energy  occurs  close  to  the  centre  of  the  vortex,  while  all  the  simulated  results  showed  a  different 
trend,  with  the  kinetic  energy  increasing  continuously  as  the  distance  to  the  wing  decreased,  with 
the  maximum  value  occuring  not  at  the  vortex  centre  but  at  the  closest  point  to  the  wing  surface. 
The  V2F  model  however  came  closest  to  predicting  the  correct  behaviour  close  to  the  wall.  This  is 
similar  to  the  result  found  by  Parneix  et  al.  [11]. 

Apsley  and  Leschziner  [10]  noted  that  RSM  offered  predictive  advantages  over  other  models 
examined  in  their  study,  especially  in  terms  of  the  far-field  structure  of  the  horseshoe  vortex.  Chen 
[13]  also  found  that  the  RSM  gave  a  superior  performance  over  that  of  the  simpler  isotropic  eddy 
viscosity  model.  In  contrast  to  this,  Rodi  et  al.  [12]  found  that  the  RSM  and  intermediate  algebraic- 
stress  and  non-linear  k-z  models  were  no  better  than  the  simpler  linear  eddy  viscosity  models  for 
the  practically  relevant  mean  quantities.  Our  own  conclusion  is  that  the  RSM  in  Fluent  can  be 
particularly  unstable  and  that  convergence  can  be  difficult  to  achieve.  Despite  this,  the  results 
shown  in  Section  7.2  calculated  using  non-equilibrium  wall  functions  on  the  coarse  grid,  even 
though  not  fully  converged,  show  agreement  with  experiment  as  good  as  that  obtained  using  the 
V2F,  Spalart-Allmaras  or  k-oo  model.  The  RSM  is  computationally  more  expensive  than  any  of  the 
other  models  however.  Both  the  Spalart-Allmaras  and  k-oo  models  require  fine  resolution  close  to 
any  boundaries,  and  this  can  also  lead  to  computationally  expensive  simulations  compared  to 
either  the  k-z  model  or  the  V2F  model.  Based  on  the  above  considerations,  it  would  appear  that  the 
V2F  model  offers  the  best  combination  of  computational  accuracy,  computational  efficiency,  and 
ease  of  use. 
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